AFIT/GA/ENY/89D-3 


STABLE  ORBITS 
ABOUT  THE  MARTIAN  MOONS 


THESIS 


Presented  to  the  Faculty  of  the  School  of  Engineering 
of  the  Air  Force  Institute  of  Technology 
Air  University 

In  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science  in  Astronautical  Engineering 


Scott  W.  Jansson,  B.S. 
Captain,  USAF 


December  1989 


Approved  for  public  release;  distribution 


unlimited . 


I 


A cknow ledge men  t  s 


I  would  like  to  thank  my  advisor,  Dr.  William  Wiesel 
for  his  expert  assistance.  His  enthusiasm  and  support  made 
this  project  both  interesting  and  educational. 


Table  of  Contents 


Page 


Acknowledgements  .  ii 

List  of  Figures  .  v 

List  of  Tables  . vii 

Not  ation  . viii 

Abstract  . .  x 

I.  Introduction  .  1 

II.  Problem  Dynamics  .  2 

Equations  of  Motion  .  4 

Kinetic  Energy  .  4 

Potential  Energy  .  6 

Potential  Energy  Due  to  Mars  ....  7 

Potential  Energy  Due  to  Moon  ....  7 

Lagrangian  .  13 

Generalized  Momenta  .  13 

Generalized  Velocities  .  14 

Hami  ltonian  .  14 

Hamilton's  Equations  .  15 

X-Axis  Symmetry  .  17 

III.  Surface  of  Section  Technique  .  18 

Restricted  Three-Body  Problem  . .  18 

Integrals  of  Motion  .  19 

Surface  .  20 

Stability  of  Orbits  .  21 


IV.  Solution  Method 


Initial  Conditions  .  22 

Trajectory  Integration  .  24 

Surface  of  Section  Points  .  24 

Orbi  t  Checks  .  2  5 

Dynamics  Verification  and  Computer  Program 

Checks  .  26 

Two-Body  Problem  .  26 

Conservation  of  the  Hamiltonian  .  27 


i  i  i 


V  .  Results  and  Discussion  .  28 

Phobos  .  28 

Typical  Orbits  .  30 

Variation  of  the  Semi-Major  Axis  .  34 

Collision  Orbits  .  35 

Orbit  Resonance  . 36 

Orbit  Evolution  .  43 

Orbit  Direction  .  46 

Deimos  .  49 

Rotating  Orbits  .  51 

Collision  Orbits  .  54 

Orbit  Evolution  .  60 

VI.  Conclusions  and  Recommendations  .  65 

Appendix  A:  Problem  Parameters  .  67 

Appendix  B:  Equations  of  Motion  .  68 

Appendix  C:  Phobos  Surface  of  Section  Plots  .  69 

Appendix  D:  Deimos  Surface  of  Section  Plots  .  93 

Bibliography  .  110 


Vita 


112 


Figure 


Page 


1.  Coordinate  System  .  3 

2.  Gravitational  Potential  Due  to  Moon  .  8 

3.  Phobos  Surface  of  Section,  H  =  -6.8528  29 

4.  Stable  Orbit  About  Phobos,  H  =  -6.8528  32 

5.  Phobos  Surface  of  Section,  H  -  -  6.852813  .  37 

6.  Phobos  Collision  Orbit,  H  =  -6.852813  .  38 

7.  Phobos  Surface  of  Section,  H  -  -  6.85287  39 

8.  Phobos  Surface  of  Section,  H  =  -  6.85287  40 

9.  Stable  Orbit  About  Phobos,  H  =  -6.8527  41 

10.  Stable  Orbit  About  Phobos,  H  =  -6.8527  42 

11.  Stable  Orbit  About  Phobos,  H  =  -6.8527  44 

12.  Stable  Orbit  About  Phobos,  H  =  -6.8527  45 

13.  Phobos  Surface  of  Section,  H  =  -6.84  47 

14.  Phobos  Surface  of  Section,  H  =  -6.84  48 

15.  Deimos  Surface  of  Section,  H  =  -2.738592  52 

16.  Deimos  Surface  of  Section,  H  =  -2.738892  53 

17.  Stable  Orbit  About  Deimos,  H  =  -2.738592  55 

18.  Deimos  Surface  of  Section,  H  =  -2.738592  56 

19.  Stable  Orbit  About  Deimos,  H  -  -2.738592  57 

20.  Deimos  Surface  of  Section,  H  =  -2.738593  58 

21.  Deimos  Collision  Orbit,  H  =  -2.738592  .  59 

22.  Deimos  Surface  of  Section,  H  =  -2.738591  61 

23.  Deimos  Surface  of  Section,  H  "  -2.73859  62 

24.  Deimos  Surface  of  Section,  H  =  -2.738589  63 

25.  Deimos  Surface  of  Section,  H  =  -2.738585  64 


Appendi x 

JL 

26  . 

Phobos 

Surface 

o  f 

Section , 

H 

- 

-6.8528  . 

.  70 

27  . 

Phobos 

Surface 

o  f 

Section, 

H 

= 

-6.8527  . 

.  71 

28  . 

Phobos 

Surface 

o  f 

Section , 

H 

- 

-6.8527  . 

.  72 

29  . 

Phobos 

Surface 

o  f 

Section  , 

H 

= 

-6.8526  . 

.  73 

30  . 

Phobos 

Surface 

o  f 

Section  , 

H 

- 

-6.8526  . 

.  74 

31  . 

Phobos 

Surface 

o  f 

Section, 

H 

= 

-6.8525  . 

.  75 

32  . 

Phobos 

Surface 

o  f 

Section, 

H 

r 

-6.8525  . 

.  76 

33  . 

Phobos 

Surface 

o  f 

Section, 

H 

= 

-6.8524  . 

.  77 

34  . 

Phobos 

Surface 

o  f 

Section  , 

H 

- 

-6.8524  . 

.  78 

35  . 

Phobos 

Surface 

o  f 

Section  , 

H 

- 

-6.8524  . 

.  79 

36  . 

Phobos 

Surface 

o  f 

Section , 

H 

-6.8524  . 

.  80 

37  . 

Phobos 

Surface 

o  f 

Section , 

H 

- 

-6.8523  . 

.  8  1 

38  . 

Phobos 

Surface 

o  f 

Sect  ion , 

H 

- 

-6.8523  . 

.  82 

39  . 

Phobos 

Surface 

o  f 

Section, 

H 

-6.8522  . 

.  83 

40  . 

Phobos 

Surface 

o  f 

Section  , 

H 

- 

-6.8522  . 

.  84 

4  1  . 

Phobos 

Surface 

o  f 

Section  , 

H 

- 

-6.8521  . 

.  85 

42  . 

Phobos 

Surface 

o  f 

Section, 

H 

-6.8521  . 

.  86 

43  . 

Phobos 

Surface 

o  f 

Sect  ion , 

H 

= 

-6.852  . 

.  87 

v 


Table 


List  of  Tables 


Page 


I  . 
II  . 


Characteristics  of 


the  Stable  Orbits  About  Phobos 
the  Stable  Orbits  About  Phobos 


33 

50 


V  L  1 


No  ta t ion 

a.b.c  axis  lengths  of  moon 

dM  differential  element  of  mass 

d  position  vector  of  the  satellite  with  respect  to 

Mars 

D  position  vector  of  the  moon  with  respect  to  Mars 

G  universal  gravitational  constant 

H  Hami 1 t  on i an 

Ijcx’Iyy’Izz  moon  mass  moments  of  inertia 

IXy,Ixz,IyZ  moon  mass  products  of  inertia 

L  Lagrangian 

m , M  ma  s  s 

generalized  momenta 
Q  ^  generalized  coordinates 

r  position  vector  of  the  satellite  with  respect 

to  the  differential  element  of  moon  mass 

R  position  vector  of  the  satellite  with  respect 

to  the  moon 

T  kinetic  energy 

v  velocity 

V  potential  energy 


X 

moon 

s 

mini mum  axis 

o  f 

i  n  e 

r  t 

i  a 

y 

moon 

s 

inte rmed i a  t  e 

a  x  i 

s  o 

f 

inertia 

z 

moon 

’  s 

maximum  axis 

o  f 

i  n  e 

r  t 

l  a 

V  1  1  1 


z 


Greek  Notation 


3  position  vector  of  differential  element  of  moon 

mass  with  respect  to  the  center  of  gravity 

w  angular  velocity  -  revolut.ion  rate  of  moon 

Q  angular  velocity  -  rotation  rate  of  moon 


AFIT/GA/ENY /89D-3 


Orbits  about  the  Martian  moons,  Phobos  and  Deimos,  were 
investigated  using  the  Poincare'  surface  of  section  tech¬ 
nique.  Hamilton's  canonical  equations  were  derived  to 
describe  the  dynamics  of  the  modified  restricted  three-body 
problem  (Mars,  moon,  artificial  satellite).  The  surface  of 
section  technique  involved  the  numerical  integration  of 
several  test  orbits  with  the  same  value  for  the  Hamiltonian. 
Apoapsis  and  periapsis  points  of  the  orbits  are  plotted  in 
the  two-dimensional  configuration  space.  Stable  orbits  were 
discovered  when  the  points  formed  sets  of  closed  curves; 
chaotic  orbits  were  indicated  by  the  scattering  of  the 
points  . 
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I  .  Introduction 

Phobos  and  Deimos,  the  Martian  moons,  are  just  two  of 
the  many  natural  satellites  in  the  solar  system  which  may  be 
modeled  as  triaxial  ellipsoids.  They  are  in  nearly  circular 
orbits  about  Mars  and  rotate  about  their  major  axes  of 
inertia.  The  objective  of  this  study  is  to  determine  wheth¬ 
er  or  not  stable  orbits  can  be  found  about  either  one  of 
these  moons . 

H  £  r  i  i  1 1  o  n  '  s  canonical  equations  will  be  derived  to 
describe  the  dynamics  of  the  modified  restricted  three-body 
problem  (Mars,  moon,  artificial  satellite).  These  equations 
will  then  be  numerically  integrated  using  the  Haming  algo¬ 
rithm  to  produce  test  orbits. 

In  the  surface  of  section  method,  numerous  test  orbits 
are  integrated  with  the  same  value  of  the  system  Hamiltoni¬ 
an.  When  a  spacecraft  passes  through  the  closest  or  fur¬ 
thest  approach  to  the  moon  (periapsis  and  apoapsis,  respec¬ 
tively),  a  point  is  plotted  in  the  orbital  plane.  If  a  set 
of  closed  curves  is  formed  by  the  points,  the  orbit  is 
stable;  otherwise,  the  orbit  is  unstable. 

Plots  of  the  periapsis  and  apoapsis  points  for  various 
Hamiltonian  values  will  be  presented.  Where  appropriate, 
stable  orbits  will  be  identified  from  these  plots. 
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I I .  Problem  Dynamics 

The  coordinate  system  used  in  this  investigation  is 
illustrated  in  Figure  1.  Phobos  and  Deimos  may  be  modeled 
as  triaxial  ellipsoids  as  shown  in  the  figure.  Each  moon 
rotates  synchronously.  That  is,  the  rotation  rate  (2)  of 
the  moon  about  its  own  axis  (z)  is  equal  to  the  revolution 
rate  (w)  of  the  moon  about  Mars.  Therefore,  the  moon's 
minimum  axis  of  inertia,  x,  remains  pointed  toward  Mars 
while  it  rotates  about  its  maximum  axis  of  inertia,  z.  The 
intermediate  axis  of  inertia,  y,  completes  the  "right-hand¬ 
ed"  set  of  orthogonal  axes. 

Mars  and  the  artificial  satellite  are  modeled  as  sym¬ 
metric  spheres  with  their  masses  concentrated  at  their 
centers.  The  center  of  mass  of  Mars  is  assumed  to  be  fixed 
in  inertial  space.  The  only  forces  assumed  to  be  acting  on 
the  three-body  system  are  the  gravitational  ones  between  the 
bodies;  all  others  are  neglected.  Additionally,  the  satel¬ 
lite  is  assumed  to  be  sufficiently  small  so  as  not  to  affect 
the  motion  of  the  other  bodies. 

The  moon's  orbit  about  Mars  is  taken  to  be  circular. 
This  assumption  is  acceptable  since  the  eccentricities  of 
Phobos  and  Deimos  are  .015  and  .00052,  respectively  (3:422). 
The  position  of  the  moon  with  respect  to  Mars  is  denoted  by 
D.  The  position  of  the  satellite  with  respect  to  the  moon 
is  denoted  by  R  and  with  respect  t.  o  Mars  by  d.  Various 


geometric  and  kinematic  parameters  used  in  this  study  for 


Mars,  Phobos ,  and  Deimos  are  given  in  Appendix  A. 

Equations  of  Motion 

The  equations  of  motion  of  the  satellite  will  be 
derived  in  the  form  of  Hamilton's  canonical  equations.  The 
Hamiltonian  is  defined  as 

H  =  -  L  (1) 

where  and  p£  are  the  generalized  coordinates  and  general¬ 
ized  momenta,  respectively.  L  is  termed  the  Lagrangian  and 
is  given  by 


L  =  T  -  V 


(2) 


where  T  and  V  are  the  kinetic  energy  and  potential  energy, 
respectively.  The  first  step  then,  in  determining  the 
equations  of  motion  is  to  find  T  and  V  in  terms  of  the 
generalized  coordinates  (Q^'s)  and  generalized  velocities 
(Qi's). 

Kinetic  Energy.  The  kinetic  energy  of  the  satellite  is 
given  by 


T 


sat 


2 


mv 


(3) 


where  m  is  the  mass  of  the  satellite  and  v  is  the  inertial 
velocity  of  the  satellite.  In  this  study,  the  inertial 
reference  frame  is  taken  to  be  at  the  center  of  mass  of 
Mars. 
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Therefore , 


v2  -  d«d  (4) 

Throughout  this  investigation,  a  dot  (•)  over  a  vector  will 
indicate  the  time  rate  of  change  of  the  vector  with  respect 
to  the  inertial  frame  (an  inertial  derivative).  As  seen  in 
Figure  1,  d  may  be  written  as 

d  =  D  +  R  (5) 

Therefore 


d  =  D  +  R  (6) 

D  may  be  expressed  as 

•  *d  rd 

D=  —  D  =  —  D+wrlxD  (7) 

dt  dt 

where  the  superscripts  i  and  r  indicate  derivatives  taken  in 
the  inertial  and  rotating  (moon-fixed)  frames,  respectively. 
Because  D  remains  constant  in  the  rotating  frame. 


rd 

-  D  =  0 

d  t 


and  Eq  (7)  simplifies  to 

D  =  w r  1  x  D 


but 


w 


r  i 


=  52  a . 


(8) 


(9) 


(  10) 
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Therefore 


D  =  (Qaz)  x  (Dax)  =  2Day  (11) 

Now  R  is  defined  to  be 

R  =  Xax  +  Yay  +  Zaz  (12) 

R  may  then  be  expressed  as 

•  M  rd  _  • 

R=— R=— R+  wr  x  R  (13) 

dt  dt 

=  Xax  +  Yay  +  Zaz  (14) 

+  (QSZ)  x  (Xax  +  Yay  +  Zaz) 

Therefore 

R  =  (X  -  2Y)ax  +  (Y  +  2X)ay  +  Zaz  (15) 

Combining  Eqs  (6),  (11),  and  (15)  results  in 

d  =  (X  -  2 Y ) a x  +  (Y  +  2X  +  2D)ay  +  Zaz  (16) 

Combining  Eqs  (3),  (4),  and  (16)  yields 

Tsat  =  ^m[(X  “  2Y)2  ♦  (Y  +  QX  +  2D)2  +  Z2  ]  (17) 

Potential  Energy.  The  potential  energy  of  the  satel¬ 
lite  may  be  divided  into  two  parts,  one  due  to  Mars  and  one 
due  to  the  moon . 

sat  mars  moon  v  ' 
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The  potential 


I  Potential  Energy  Due  to  Mars. 

energy  of  the  satellite  due  to  the  gravitational  attraction 
of  Mars  is  easily  derived.  In  this  study.  Mars  is  assumed 
to  possess  spherical  symmetry.  Therefore,  the  center  of 
mass  of  Mars  is  assumed  to  coincide  with  its  center  of 
gravity.  The  potential  energy  of  the  satellite  due  to  the 
gravitational  attraction  of  Mars  may,  therefore,  be  written 
in  the  simple  Newtonian  form  as 


V 


ma  r  s 


GM 


ma  r  s 


m 


d 


where  d  is 


d=  | d |  =  |D+R|  =  [(X  +  D)2  +  Y2  +  Z2]1/2 


(19) 


(20) 


Potential  Energy  Due  to  Moon.  The  potential 
energy  of  the  satellite  due  to  the  gravitational  attraction 
of  the  moon  is  more  difficult  to  determine  since  the  moon  is 
not  spherically  symmetric.  Its  shape  is  assumed  to  be  a 
triaxial  ellipsoid  as  shown  in  Figure  2.  The  gravitational 
potential  for  an  arbitrarily  shaped  body  is  given  in  Meiro- 
vitch  (9:430-436). 


The  potential  energy  in  an  inverse  square  force  field 
may  be  written  as 


V 


moon 


rdM_n 

I  mo  o  n 

Gm - 


(21) 
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The  scalar  r  is  the  distance  between  the  satellite 
center  of  gravity  and  a  mass  differential  element  of  the 
moon,  dMjjoon'  It  may  be  written  as 

r  =  |  r  j  =  |R  -  (3 1  (22) 

Its  inverse  may  be  written  as 


r-1  =  |R  -  p|_1  =  [(R  -  P)2]~1/2  (23) 


(R2  - 

2R«p  +  0 2)-1/2 

2 

(24) 

[  R2  <  1 

-  2———  +  ^-)]-l/2 

R2  R2 

2 

(25) 

R  ~ 1  [  1 

-  (2^1  +  5_)]-l/2 

R2  R2 

(26) 

Using  a  binomial  expansion,  it  can  then  be  written  in  the 
form  of  a  power  series  as 


_ - 1 , ,  R*P  1  .  P  .  2  3  R»P  t 8  \  2  i  2 

r2  2  R  8  r2  R 

=  R'1  +  R'3(R«0)  -  ^R"3P2  +  -^R-5  ( R»  3 ) 2 

2  2 


-  -R~5(R*p)p2  +  -R_5p4  +  ... 
2  8 


(27) 

(28) 


From  Figure  2,  the  following  vectors  are  defined 


R  =  Xax  +  Yay  +  Zaz  (29) 

P  =  aax  +  bay  +  caz  (30) 
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Therefore 


R2  =  X2  +  Y2  +  Z2 


l2  =  a2  +  b2  +  c2 


(31) 

(32) 


R»  3  =  aX  +  bY  +  cZ 


(33) 


Combining  Eqs  (21),  (28),  and  (31)-(33)  results  in 


vmoon  =  -  O-O’1 


dMmoon  “  GmR 


-3 


(aX  +  bY  +  cZ)dMmoon  (34) 


+  —GmR 
2 


-3 


(a2  +  b2  +  c2 )dM 


moon 


3  -s 
-GmR  5 


(aX  +  bY  +  c  Z ) 2  dM 


moon 


The  higher  order  terms  resulting  from  the  binomial  expansion 
have  been  dropped  since  their  effect  is  negligible. 

Eq  (34)  may  now  be  broken  down  into  smaller  parts. 


-  GmR 


-1 


d^raoon  GmR  ^moon 


(35) 


Because  the  origin  is  taken  to  be  at  the  center  of  mass, 


-  GmR 


-3 


(aX  +  bY  ♦  cZ)dMmoon  =  0 


(36) 


Combining  the  third  and  fourth  terms  of  Eq  (34)  results  in 


1  _  3 

-GmR  3 


(3R  2 ( aX  +  bY  +  cZ)2  -  (a2  ♦  b2  +  c2)]dMmoon  (37) 


1  _  3 
-GmR  3 


[ 3R-2 ( a2X2  +  b2  Y2  +  c2Z2) 


(38) 
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+  6 ( abXY  +  acXZ  +  bcYZ)  -  (a2  +  b2  +  c2)]dM 


moon 


■jGmR"3[(3R  2X2  -  1) 


a2dM 


moon 


(39) 


( 3R~2 Y2  -  1) 


b2dMmoon  +  ( 3R-2  Z2  -  1) 


c2dM 


moon 


+  6  (  XY 


abdMmoon  +  XZ 


acdMmoon  +  YX 


bcdMmoon^ 


But 


a^dM 


moon 


[(a2  +  b2)  +  (a2  +  c2)  -  (b2  +  c2)]dM 


moon 


b2dM 


moon 


X 

"2 ' 1  z  z  +  xyy  ~  1  xx  ^ 


[(a2  +  b2)  +  (b2  +  c2)  -  (a2  +  c2 ) ] dM  (42) 


(40) 


(41) 


moon 


X 

2^  *zz  +  *xx  ~  Jyy ) 


(43) 


c2dM 


moon 


((a2  +  c2)  +  (b2  +  c2)  -  (a2  +  b2)]dM 


moon 


(44) 


2( xyy  +  Ixx  “  1  z  z  ^ 


(45) 


abdM  =  I 

moon  Axy 


(46) 


acdMmoon  *xz 


(47) 


bcdr^moon  ^  y  z 


(48) 


1 1 


Ixx,  I  y  y ,  and  Izz  are  the  mass  moments  of  inertia.  I  xy 
Ixz,  and  IyZ  are  the  mass  products  of  inertia.  Using  prin¬ 
cipal  axes,  the  mass  products  of  inertia  for  the  triaxial 
ellipsoid  are 


*xy  *xz  ~  Iyz  0 


(49) 


The  mass  moments  of  inertia  are  calculated  using  the  follow¬ 
ing  equations  (8:542) 


xx 


5«»oon<b^  * 


'yy 


5^moon^a  +  c  ) 


z  z 


j^moon^3  +  b  ) 


where  a,  b,  and  c  are  the  axis  lengths  of  the  moon. 

The  potential  energy  due  to  the  moon  may  now  be  written 


moon 


-  GmMmoonR 


-1 


-  ^GmR‘3[ (3R~2X2  -  l)(lyy  +  Izz 


*xx  ) 


(  50) 


+  ( 3R  Y  -  1)(IXX  +  Izz  -  Iyy) 


+  ( 3R_2Z 2  -  1)(IXX  +  Ivv  -  IZ2) 1 


moon 


=  -  GmM 


moon 


R-'  *  iG„R-3(Ixx  ♦  Iyy  +  I„> 


-  — GmR-  3  [  X3  (  Iyy  ♦  I„  -  Ixx) 


(51) 


*  1  *  1zz  '  Iv»>  -  ZZU„  *  I.v  -  I,,)] 


yy 


xx  *  y  y 


z  z 
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The  next  step  in  determining  the  equations 


of  motion  is  to  find  the  Lagrangian,  Combining  Eqs  (2), 
(17),  (19),  (20)  and  (51)  the  Lagrangian  becomes 

L  =  T  -  V  =  ^m[(X  -  2  Y ) 2  +  (Y  +  2X  +  2D)2  +  Z2]  (52) 

*  “aoo."11'1  -  XG”R~3(I**  *  lyy  *  Izi> 

*  “GmR  [X  (lyy  +  ^22  ~  ^ 

*  l2(I„  *  Izz  -  lyy)  *  z2(Ixx  *  lyy  -  ‘nil 

Dividing  out  the  mass  of  the  satellite,  m,  the  Lagrangian 
may  be  written  on  a  per  unit  mass  basis  as 

L  =  ^((X  -  QY)2  +  (Y  +  2X  +  2D)2  +  Z2 ]  (53) 

+  Omars'*'1 

*  '»..onR'“  -  XCR'3(Ix*  ♦  lyy  +  G*) 


+  -GR"5(X2 
4 

^  Iyy  + 

*  z  z 

Ixx> 

+  Y2 (  I  + 

^  z  z 

iyy) 

♦  z2<Iyx  -  1 

yy  ~  ^zz^l 

Generalized  Momenta. 

Next, 

the 

gene  ra 1 i zed 

mome n t  a  , 

arp  derived 

PX  = 

6  L 

— r  =  X 
oX 

-  2Y 

(  54) 

PY  = 

oL 

— r  =  Y 
6  Y 

+  2X 

+  2D 

(55) 
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p 


z 


6  L 
oZ 


(56) 


Generalized  Velocities.  Before  the  Hamiltonian  can  be 
defined,  the  generalized  velocities,  's,  must  first  be 
eliminated  from  the  Lagrangian.  Rearranging  Eqs  (54)-(56) 
results  in 


X  = 

px  + 

2Y 

(57) 

Y  = 

P  Y  - 

2X  - 

2D 

(58) 

Z  = 

Pz 

(59) 

Substituting  Eqs 

(54)— (56) 

into 

the 

Lagrangian 

(Eq 

(53)) 

provides 

1 

L  =  - 
2 

(px2  *  py2 

+  pz2 

) 

(60) 

+ 

+ 

1 

-  -GR 
4 

"3(  I 

+  1  4* 

xx  Lyy 

+ 

-gr"5[x2(i 

4 

+ 

yy 

*  z  z 

-  rx«) 

+ 

*  I 

z  z 

1  y  y ) 

*  z2('xx  * 

lyy 

-  Izz>l 

Hami ltonian  .  With  the  Lagrangian  in  the  correct  form, 
the  Hamiltonian  may  be  written  as  follows 

H  =  -  L  =  PX(PX  +  2Y)  +  Py(PY  “  2X  -  2D)  +  Pz2  (61) 

"  ^(pX2  +  PY2  +  pZ2) 
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-  GMm_  _cd” 1 
ma  r  s 


-  '".o.n*'1  *  ^CR"3<-Xx  *  !yy  *  I2z> 

-  jOR-5[X2(Iyy  *  I„  -  I„) 


*  »2»„  *  !Z2  -  lyy)  *  !2(I„  *  Ivv  -  I*,)] 


xx  yy  z  z 


-(Px2  +  Py2  +  Pz2)  +  PXQY  -  Py2(X  +  D)  (62) 


-  GMmarsd' 


-  GMmnnnR_1  +  T<*  3Uvv  +  IVv  +  I,,) 


xx  yy  z  z 


4GR  ^yy  +  Xzz  *xx^ 


+  y2(rxx  +  lzz  -  lyy)  +  z2(Ixx  +  ryy  "  lzz^ 


Hamilton's  Equations.  The  equations  of  motion  may  now 
be  written  in  canonical  form.  It  is  advantageous  to  use 
this  form  because  the  resulting  first  order  differential 
equations  can  be  numerically  integrated  more  readily.  In 
addition,  since  the  Hamiltonian  contains  only  the  general¬ 
ized  coordinates  and  momenta,  and  not  their  time  deriva¬ 
tives,  first  integrals  of  motion  are  easier  to  determine. 


(63) 


=  Pv  +  2Y 


(  64) 


Y  =  -  =  Pv  -  2(X  +  D) 

oP  v  T 


(65) 


oH 

Z  =  —  = 

oP, 


oH 

oQl 


P 


X 


oH 

oX 


PY 


oH 

6Y 


These  equations 
convenience . 


(66) 


(67) 


=  2PV  -  GMmarsd"3(X  +  D)  -  GMm„„nR‘3X  (68) 


moon 


♦  -ClT5X(3Iyy  *  3I„  -  Ixx) 


-  ^GR-7XtX2(Iyy  +  I„  -  Ixx) 


*  ^(I„  -  Iz2  -  :yy)  *  ZMIxx  *  !yy  - 


=  -  2PX  -  GMmarsd-2Y  -  GMboooR*3Y 


(70) 


3  _c  , 

+  -GR  5 Y ( 31 Y  +  31  -  Ivv) 


XX 


zz  yy 


-  i^GR-7Y[X2(Iyy  *  I„  -  Ixx) 


+  Y*'(I,  M  T  -  T  '  <*■  *2 , 


xx 


i7z  -  iyy )  ♦  z^(ixx  ♦  I 


yy 


=  -  GMmarsd'3z  “  GMmoonR"3z 


(70) 


+  4GR_5z(3Ixx  +  3Iyy  "  I2z> 


4  GR  Z[X  (Iyy  +  Izz  "  Ixx) 


+  ^  ^xx  +  ^zz  -  ^yy^  +  Z*”^xx  +  *yy 


Izz)l 


of  motion  are  also  listed  in  Appendix  B  for 
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X-Axis  Symmetry 

The  equations  of  motion  may  be  transformed  using  the 
following  relationships 


X - >  X 

Y - >  -Y 

t - >  -t 

When  this  transformation  is  made,  »-he  equations  of  motion  do 
not  change  in  the  rotating  coordinate  system.  Therefore, 
for  every  orbit,  there  is  a  second  orbit  which  is  a  reflec¬ 
tion  of  the  first  orbit  about  the  XZ  plane.  The  second 
orbit,  however,  is  traversed  in  the  opposite  direction  from 
the  first  due  to  -t .  This  symmetry  becomes  very  useful  in 
the  surface  of  section  technique  and  effectively  doubles  the 
number  of  data  points  obtained  for  a  given  trajectory. 
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Poincare's  surface  of  section  method  is  used  to  discov¬ 


er  stable  orbits  around  the  Martian  moons.  In  this 
technique,  the  dynamics  of  the  three-body  problem  are  re¬ 
stricted  to  the  planar  case  and  the  one  known  integral  of 
motion  is  determined.  A  surface  of  section  is  then  intro¬ 
duced  into  the  problem  in  order  to  determine  if  further 
integrals  of  motion  exist.  If  a  second  integral  of  motion 
exists,  stable  orbits  may  be  discovered. 

Restricted  Three-Body  Problem 

In  what  is  commonly  referred  to  as  the  Restricted 
Three-body  Problem,  only  the  gravitational  forces  between 
the  bodies  are  considered.  Perturbations  due  to  other 
forces  such  as  the  gravitational  attraction  of  other  bodies, 
atmospheric  drag,  a  non-spher ical  body,  etc.  are  assumed  to 
be  very  small  and  are  ignored.  The  equations  previously 
derived  describe  the  motion  of  a  satellite  in  a  Modified 
Restricted  Three-body  Problem.  Gravity  terms  due  to  the 
non-s phe r i ca 1  moon  are  also  included. 

The  artificial  satellite  in  the  restricted  three-body 
problem  is  constrained  to  motion  in  the  plane  formed  by  the 
moon  revolving  around  the  planet.  Because  of  this  planar 
restriction,  Eqs  (66)  and  (70)  can  be  ignored  and  Z  can  be 
set  equal  to  0  in  the  remaining  equations  (Eqs  (64),  (65), 

(68),  and  (69)).  These  equations  of  motion  describe  a 


system  with  two  degrees  of  freedom  and  four  dimensions  in 


the  phase  space  . 


Integrals  of  Motion 

When  four  initial  conditions  are  defined,  a  solution  of 
the  system  may  be  represented  in  the  four-dimensional  phase 
space.  If  an  integral  of  motion  exists,  solutions  of  the 
system  may  be  represented  in  a  three-dimensional  subspace 
for  particular  values  of  the  constant  integral  of  motion 
(13:127). 

In  this  problem,  the  Lagrangian  (Eq  60)  does  not  depend 
explicitly  on  time.  Therefore 


oL 

—  =  0 
6 1 


(71) 


An  integral  of  motion,  known  as  Jacobi's  integral  may  then 
be  determined  from  the  following  (9:83) 


E 


L  -  constant 


(72) 


In  this  case,  the  Jacobi  integral  turns  out  to  be  equal  to 
the  system  Hamiltonian  (Eq  62).  This  was  to  be  expected 
since  the  Hamiltonian  does  not  depend  explicitly  on  time. 

In  this  problem,  the  artificial  satellite  is  con¬ 
strained  to  move  on  a  three-dimensional  manifold  within  the 
four-dimensional  phase  space  because  of  the  existence  of  the 
Jacobi  Integral  (7:6).  The  Hamiltonian  provides  a  one¬ 
dimensional  pa r ame t r i z a t i on  of  the  manifolds  on  which  the 


satellite  motion  takes  place. 

If  another  independent  integral  of  motion  exists,  the 
satellite  is  further  constrained  to  move  on  a  two- 
dimensional  manifold  embedded  in  the  three-dimensional  one 
defined  by  the  Hamiltonian  (7:6).  The  surface  of  section 
technique  is  used  to  determine  the  existence  of  such  inde¬ 
pendent  integrals  of  motion. 

Surf  ace 

In  the  surface  of  section  method,  a  surface  is  intro¬ 
duced  into  the  phase  space  according  to  some  given  relation¬ 
ship.  If  a  second  integral  of  motion  exists,  the  intersec¬ 
tion  of  the  two-dimensional  manifold  on  which  the  satellite 
is  constrained  to  move  and  the  introduced  surface  will 
generally  be  one-dimensional  (7:6).  If  a  second  integral  of 
motion  does  not  exist,  the  intersection  of  the  surface  and 
the  satellite  orbit  will  not  be  one -d imens i ona 1  .  Instead, 
the  intersection  will  be  scattered  over  a  larger  subset  of 
the  surface  . 

In  this  study,  the  surface  of  section  is  defined  by  the 
following  condition 

R  •  V  =  0  (73) 

or  for  the  planar  case  under  consideration 

XX  +  YY  =  0  (74) 

This  condition  occurs  at  the  satellite  orbit's  apoapsis  and 
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periapsis  points.  These  are  the  points  where  the  satellite 
is  farthest  from  and  nearest  to  the  moon.  Points  along  the 
orbit  of  the  satellite  where  this  condition  occurs  are 
plotted  in  the  two-dimensional  configuration  space. 

Other  criteria  for  the  surface  of  section  are  possible; 
however,  R»V  =  0  is  used  because  it  indicates,  by  inspec¬ 
tion,  whether  or  not  an  orbit  will  collide  with  the  moon. 
If  any  of  the  surface  of  section  points  lie  beneath  its 
surface,  the  trajectory  will  obviously  collide  with  the 
moon  . 

Stability  of  Orbits 

If  the  surface  of  section  points  plotted  in  the  config- 

* 

uration  space  form  closed  contours,  a  stable  orbit  is 
present  for  that  particular  value  of  the  Hamiltonian.  A 
periodic  orbit  will  pass  through  the  centers  of  the  closed 
curves  . 

If  the  surface  of  section  points  form  closed  curves 
that  intersect  one  another,  unstable  orbits  may  be  located 
which  pass  through  the  intersection  points.  The  intersec¬ 
tion  of  the  curves  is  often  referred  to  as  a  saddle  point 
(an  unstable  equilibrium  point). 
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IV  .  Solution  Method 

Satellite  trajectories  are  determined  through  numerical 
integration  with  a  FORTRAN  computer  program.  The  computer 
program  allows  for  significant  interaction  between  the 
operator  and  computer. 

The  initial  X  and  Y  coordinates  for  the  trajectory  and 
the  Hamiltonian  value  are  the  primary  inputs  along  with  the 
integration  step  size  and  the  total  integration  time.  A 
data  file  with  the  X  and  Y  coordinates  for  all  of  the  sur¬ 
face  of  section  points  is  the  main  output. 

The  output  data  files  for  several  trajectories,  each 
having  the  same  Hamiltonian  value,  are  then  plotted  with  a 
two-dimensional  graphics  package.  From  these  plots,  the 
discovery  of  orbits,  both  stable  and  unstable,  may  be  made. 
Initial  Conditions 

Prior  to  the  integration  of  the  equations  of  motion, 
initial  conditions  for  the  states  X,  Y,  P and  P  y  must  be 
determined.  The  initial  values  of  X,  Y,  and  the  Hamiltoni¬ 
an,  H,  are  chosen.  In  addition,  the  starting  point  of  all 
trajectories  is  chosen  so  the  surface  of  section  condition 
of  Eqs  (73)  and  (74)  is  met.  Given  these  conditions,  the 
initial  values  for  P^  and  Py  may  be  calculated  by  first 
obtaining  a  quadratic  expression  for  Y. 

To  obtain  a  quadratic  expression  for  Y,  the  following 
equations  are  utilized 
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oL 

PX  =  ^  =  x  -  =v 
6L 

Pv  =  — r  =  Y  +  2X  +  2D 
Y  oY 


(54) 

(55) 


1 

2 


H  =  -(Pv2  +  Pv2  +  P?2)  +  Py2Y  ~  PV2(X  +  D)  (62) 


-  GMmarsd 


-1 


GMmoon^  +  4G^  ^xx  +  *yy  +  *zz^ 
— (JR  ^[X2(Iyy  +  I2Z  ”  Ixx) 


+  Y2(IXX  +  IZZ  -  Iyy)  +  Z2(IXV  +  Ivv  -  177  )] 


xx  y y  z  z 


XX  +  YY  =  0 


(74) 


Substituting  Eqs  (54)  and  (55)  into  Eq  (62),  and  combining 
all  of  the  potential  energy  terms  into  one  quantity  (V), 
results  in  the  following  expression  for  the  Hamiltonian 

H  =  Ux2  +  Y2  -  22X2  -  Q2Y2  -  22  D2  )  (75) 

2 

-22YX  -  Q2DX  +  V 
Solving  Eq  (74)  for  X  provides 

Y  • 

X  =  -  -Y  (76) 

X 

Substituting  Eq  (76)  into  Eq  (75)  yields  the  following 
quadratic  expression  for  Y 
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(77) 


(X2  +  Y2  +  4QY2 ) Y2  +  ( 42XY  2 ) Y 

+  (2X2)["(Q2X2  +  22 Y2  +  Q2D2)  -  22DX  +  V  -  H]Y  =  0 

Solving  the  above  expression  for  Y  results,  of  course, 
in  two  different  values.  Either  value  may  be  chosen.  In 
many  instances,  choosing  one  value  will  result  in  a  prograde 
trajectory,  while  choosing  the  other  value  will  result  in  a 
retrograde  trajectory. 

X  may  be  calculated  from  Eq  (76).  The  required  initial 
values  for  the  generalized  momenta,  and  Py,  may  then  be 

determined  from  Eqs  (54)  and  (55). 

Trajectory  Integration 

Haming's  Ordinary  Differential  Equations  Integrator 
is  used  to  integrate  the  equations  of  motion  (Eqs  (64), 
(65),  (68),  and  (69)).  It  is  a  fourth  order  predictor- 

corrector  algorithm.  It  extrapolates  the  last  four  values 
of  the  state  vector  to  obtain  a  predicted  next  value  (the 
prediction  step)  (16:120).  It  then  evaluates  the  equations 
of  motion  at  the  predicted  value  and  corrects  the  extrapo¬ 
lated  point  using  a  higher  order  polynomial  (the  correction 
step)  (16:120).  The  Haming  method  is  numerically  very 
stable  and  provides  high  accuracy. 

Surface  of  Section  Points 

After  each  integration  step,  the  value  of  R*V  is  calcu¬ 
lated.  The  current  value  is  compared  to  the  last  value  and 


if  the  sign  of  the  value  has  changed,  the  surface  of  section 
condition  (Eq  73)  has  been  met. 

In  order  to  obtain  very  precise  values  for  the  surface 
of  section  points,  the  last  four  values  of  the  integration 
time,  the  X  and  Y  coordinates,  and  the  R»V  calculation  are 
maintained.  When  a  check  of  the  R»V  value  indicates  the 
orbit  has  passed  through  a  surface  of  section  point,  the 
Newton  interpolation  polynomial  (1:112)  is  used  with  the 
last  four  integration  states  to  obtain  a  very  precise  deter¬ 
mination  of  the  coordinates  of  the  surface  of  section  point. 
R*V  is  a  function  of  X  and  Y,  however.  Therefore,  inverse 
interpolation  is  used  to  determine  the  X  and  Y  coordinates 
for  which  R  •  V  =  0  (1:119). 

As  discussed  earlier,  there  is  trajectory  symmetry 
about  the  X-axis.  Therefore,  in  addition  to  the  surface  of 
section  point  determined  from  the  trajectory  being  inte¬ 
grated  ( X , Y ) ,  the  reflection  of  that  point  about  the  X-axis 
( X , -Y )  is  also  obtained. 

Orbit  Checks 

Several  checks  are  made  during  the  integration  of  the 
trajectories.  The  location  of  the  satellite  relative  to  the 
moon  is  continuously  checked. 

If  the  satellite  travels  too  far  away  from  the  moon  an 
escape  trajectory  is  assumed.  In  many  cases,  the  satellite 
escapes  into  a  Martian  orbit.  If,  on  the  other  hand,  the 
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satellite  passes  too  closely  to  the  center  of  the  moon,  the 
singularity  at  that  point  causes  the  integration  to  "blow 
up."  Whenever  either  one  of  these  conditions  occurs,  inte¬ 
gration  of  the  trajectory  is  halted. 

Even  though  trajectories  which  pass  through  (or  very 
close  to)  the  center  of  the  moon  become  numerically  unsta¬ 
ble,  other  trajectories  continue  right  through  collisions 
with  the  moon  without  difficulty.  These  trajectories  are 
allowed  to  continue,  since  they  may  provide  useful  results 
in  the  surface  of  section  technique.  In  any  event,  a  check 
of  the  trajectory  will  determine  that  a  "collision"  has 
occurred  in  these  cases  and  make  that  fact  known. 

Dynamics  Verif ' _ ^.ion  and  Computer  Program  Checks 

In  order  to  verify  that  the  equations  of  motion  were 
correctly  derived  and  that  the  computer  program  was  perform¬ 
ing  as  intended,  various  checks  were  made. 

Two-Body  Problem.  One  verification  check  involved 
reducing  the  dynamics  of  the  problem  to  a  simple  two-body 
problem.  In  this  case,  the  gravitational  attraction  of  Mars 
(GMmars),  the  system  rotation  (2),  the  moon's  moments  of 
inertia  (Ixx.  ^yy*  anc*  -  and  the  distance  to  the  origin 

of  the  coordinate  system  (D)  were  all  set  to  0.  The  satel¬ 
lite  was  then  given  circular  speed  relative  to  the  moon  for 
its  initial  altitude.  The  orbit  integrated  was  indeed 
circular  and  did  return  to  the  original  starting  point. 
When  the  moments  of  inertia  were  included,  the  orbit  became 
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slightly  elliptical  and  contracted  as  expected. 

In  another  case,  the  gravitational  attraction  of  the 
moon  (GMmoon)  anc*  the  initial  velocity  of  the  satellite  were 
set  to  0.  The  starting  point  of  the  trajectory  was  chosen 
very  near  to  the  moon.  The  satellite  remained  relatively 
stationary  in  the  rotating  coordinate  system,  indicating  it 
was  indeed  in  the  expected  circular  orbit  of  mars. 

Conservation  of  the  Hamiltonian.  Another  check  in¬ 
volves  the  Hamiltonian.  As  discussed  previously,  the  Hamil¬ 
tonian  is  a  constant  integral  of  motion.  It  should  remain 
invariant  over  the  entire  trajectory.  H,  therefore,  is 
periodically  checked  throughout  the  trajectory  integration 
to  insure  it  is  conserved  to  several  decimal  places. 
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V •  Results  and  Discussion 

Stable  orbits  were  discovered  about  each  of  the  Martian 
moons.  The  discussion  of  results  will  be  presented  sepa¬ 
rately  for  Phobos  and  Deiraos  .  However,  the  nature  of  the 
orbits  about  each  of  the  moons  is  quite  similar,  as  might  be 
expec  ted . 

Phobos 

Figure  3  is  a  typical  surface  of  section  plot  for 
Phobos.  Each  of  the  test  trajectories  was  initiated  along 
the  negative  X-axis  (between  the  moon  and  Mars).  Several 
stable  trajectories  (3-9)  were  integrated  for  each  chosen 
Hamiltonian  value.  They  were  evenly  spaced  out  over  the 
"region  of  stability".  In  this  region,  all  of  the  test 
trajectories  remain  in  orbit  around  Phobos  for  the  entire 
integration  time.  They  demonstrate,  at  the  very  least, 
practical  stability.  The  trajectories  initiated  beyond  this 
region  were  unstable.  Their  surface  of  section  points  were 
scattered  in  the  configuration  space. 

In  Figure  3,  eight  stable  test  trajectories  were  inte¬ 
grated  for  the  Hamiltonian  value,  H  =  -6.8528.  Separate 

curves  for  each  of  them  are  readily  identifiable  in  the 
figure. 

Trajectories  which  were  initiated  outside  the  stable 
region,  closer  to  the  moon,  proceeded  to  the  moon's  center. 
As  expected,  unstable  trajectories  closer  to  the  moon 
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reached  the  center  in  far  fewer  orbits  than  those  further 


away  from  the  moon.  Trajectories  that  were  initiated  out¬ 
side  the  stable  region,  further  from  the  moon,  escaped  the 
moon's  influence.  They  normally  entered  a  Mars  orbit  simi¬ 
lar  to  the  orbit  of  the  moon  itself.  These  characteristics 
were  similar  to  those  of  previous  thesis  efforts. 

All  of  the  trajectories  were  integrated  for  5,000,000 
seconds  (57.9  days).  This  integration  time  provided  thou¬ 
sands  of  surface  of  section  points!  enough  so  the  size  and 
shape  of  the  curves  formed  were  easily  recognizable.  The 
integration  step  size  normally  used  was  200  seconds.  Under 
these  conditions,  the  run  time  on  a  main  frame  computer  was 
several  minutes.  In  a  few  instances,  a  smaller  step  size 
(50-100  seconds)  was  required  to  prevent  numerical  integra¬ 
tion  problems!  usually  for  trajectories  close  to  the  center 
of  the  moon.  Even  with  the  large  step  size,  the  location  of 
the  surface  of  section  points  was  accurately  calculated  with 
the  inverse  interpolation  method  discussed  earlier. 

Typical  Orbits.  For  the  typical  orbit,  four  points 
meeting  the  surface  of  section  criteria  (R»V  =  0)  were 

determined  each  time  a  trajectory  encircled  the  moon.  For 
the  retrograde  trajectories  represented,  the  points  appeared 
successively  in  a  counter-clockwise  manner.  The  four  re¬ 
gions  where  the  points  were  determined  are  easily  identifia¬ 
ble  in  Figure  3.  The  curves  formed  for  each  trajectory  are 
closed  and  thus  indicate  stable  orbits.  However,  most  of 
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the  orbits  rotate  or  precess  about  the  moon  and  are  not 
simple,  closed  periodic  orbits.  The  four  surface  of  section 
points  rotate  regularly  within  their  "own  region"  for  each 
successive  encirclement  of  the  moon. 

The  centers  of  the  closed  curves  represent  the  surface 
of  section  points  for  a  stable,  closed  periodic  orbit. 
Figure  4  shows  the  simple  periodic  orbit  associated  with  the 
surface  of  section  of  Figure  3.  In  order  to  demonstrate 
that  the  surface  of  section  points  for  the  closed  periodic 
orbit  do  in  fact  coincide  with  the  centers  of  the  closed 
curves,  the  orbit  has  been  superimposed  with  the  surface  of 
section  plot.  The  orbit's  eccentricity  is  .221  and  its 
period  is  13,900  seconds  (3.86  hours).  A  satellite  placed 
in  this  orbit  would  encircle  the  moon  approximately  twice 
during  the  orbital  period  of  Phobos  about  Mars. 

The  surface  of  section  technique  was  used  to  discover 
many  stable  orbits  about  Phobos.  The  orbit  altitudes  vary 
from  just  above  the  moon's  surface  to  several  hundred  kilom¬ 
eters  away.  The  surface  of  section  plot  associated  with 
each  of  the  orbits  is  contained  in  Appendix  C.  They  are 
each  identified  with  their  particular  Hamiltonian  value,  H. 

Table  I  provides  a  summary  of  the  characteristics  of 
each  of  the  stable,  closed  periodic  orbits.  The  period  of 
the  first  orbit  listed  is  approximately  half  that  of  Pho- 
bos's  orbit  about  Mars.  The  period  of  all  of  the  other 
orbits  is  approximately  equal  to  Phobos ' s  orbit  period 
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Figure  4.  Stable  Orbit  About  Phobos,  H  z  -6.8528 
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Table  I.  Characteristics  of  the  Stable  Orbits  About  Phobos 
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86 
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132 
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217 
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(Phobos  rotates  beneath  each  orbit  in  the  rotating  coordi¬ 
nate  system).  The  eccentricities  of  the  orbits  are  all 
about  .34  -  .35  with  the  exception  of  the  first  orbit  which 

is  slightly  less  eccentric  (.221).  The  table  also  indicates 
the  altitudes  of  the  orbits  increase  at  a  lower  and  lower 
rate  as  the  value  of  H  increases. 

The  initial  velocities  of  the  orbits  range  from  11.4 
m/s  for  the  lowest  orbit  (H  =  -6.8528  )  to  327  m/s  for  the 
highest  orbit  (H  =  -6.84). 

Variation  of  the  Semi-Major  Axis.  Figures  3  and  4 
demonstrate  an  interesting  and  unexpected  phenomenon.  For 
each  orbit  in  a  normal  Keplerian  ellipse,  only  two  points 
will  meet  the  surface  of  section  criteria  (R»V  -  0)  used  in 
this  investigation.  These  two  points,  apoapsis  and  periap- 
sis,  are  found  at  the  ends  of  the  major  axis.  For  the 
closed  orbit  of  Figures  3  and  4,  however,  four  points  meet¬ 
ing  the  surface  of  section  criteria  were  discovered.  This 
orbit  does  not  represent  a  simple  rotating  Keplerian  ellipse 
as  expected. 

The  semi-major  axis  is  not  constant  as  the  satellite 
orbits  the  moon.  Instead,  it  varies  with  a  period  equal  to 
half  that  of  the  satellite's  orbital  period.  It.  can  be  seen 
from  Figure  4  that  a  second,  longer  apoapsis  distance  is 
achieved  during  the  first  quarter  of  the  orbit  period. 
However,  after  half  of  the  orbit  period,  the  major  axis  has 
shifted  back  to  its  original  position  and  the  expected 
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periapsis  point  is  achieved. 

During  the  second  half  of  the  orbit,  the  major  axis 
shifts  again.  The  resulting  apoaps is  is  symmetric  about  the 
X-axis  with  respect  to  the  larger  apoapsis  which  occurred 
during  the  first  half  of  the  orbit  period.  By  the  time  the 
satellite  has  completed  one  orbit,  the  major  axis  has  shift¬ 
ed  back  to  its  original  position  again.  This  perturbation 
occurs  in  all  of  the  orbits  discovered  around  Phobos . 

This  phenomenon  was  discovered  by  Tycho  Brahe  in  the 
orbit  of  the  Earth's  moon  (  1  1  :  289  ).  Brahe  found  that  the 
semi-major  axis  of  the  Moon's  orbit  executed  small  oscilla¬ 
tions.  This  perturbation,  termed  ’’the  variation,"  has  a 
period  of  half  a  synodic  month.  It  results  from  the  fact 
that  there  is  less  pull  on  the  Moon  in  the  radial  direction 
of  the  Earth  at  new  and  full  moon  (the  force  potential  is  a 
maxima)  than  at  the  quarters  (the  force  potential  is  a 
minima)  (2:287).  The  curvature  of  the  Moons ’ s  orbit,  there¬ 
fore,  is  least  at  new  and  full  moon  and  greatest  at  the 
quarters,  so  that  the  orbit  is  elongated,  with  its  longer 
axis  perpendicular  to  the  direction  of  the  Sun  (2:287). 

The  results  obtained  in  this  study  demonstrate  the  same 
perturbation  phenomenon.  In  this  case.  Mars  is  the  third 
body  that  causes  the  effect.  As  in  the  case  of  the  Moon's 
orbit,  the  satellite  s  major  axis  in  this  study  is  elongated 
in  a  direction  perpendicular  to  the  third  body. 

Collision  Orbits.  As  discussed  earlier,  the  surface  of 
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section  criteria  provides  a  clear  indication  of  orbit  tra¬ 
jectories  that  will  impact  the  moon.  Figure  5  clearly 
indicates  all  of  the  trajectories  integrated  are  collision 
orbits.  All  of  the  periapsis  points,  and  in  some  cases  even 
the  apoapsis  points,  lie  below  the  surface  of  the  moon.  The 
stable,  closed  periodic  orbit  indicated  does  in  fact  lie 
partially  below  the  moon's  surface  as  seen  in  Figure  6. 

Orbit  Resonance.  Nonlinear  systems  often  exhibit 
resonances  at  any  rational  multiple  of  the  forcing  frequency 
(16:140).  In  the  case  at  hand,  resonances  occur  when  the 
period  of  the  satellite  is  a  rational  multiple  of  the  moon's 
orbit  period.  Resonant  orbits  usually  produce  several 
"islands"  associated  with  a  stable  periodic  orbit.  These 
"chains  of  islands"  were  first  described  by  Henon  and  Heiles 
(5:76)  . 

Figures  7  and  8  provide  an  excellent  indication  of 
orbit  resonance.  Figure  7  shows  the  entire  surface  of 
section  plot  for  H  =  -6.8527.  Figure  8  is  a  magnification 
of  the  same  surface  of  section  plot  for  the  region  closest 
t o  Mar s . 

As  expected,  a  stable,  closed  periodic  orbit  is  located 
by  the  center  of  the  four  concentric  closed  curves  of  Figure 
8.  Figures  9  and  10  show  the  simple  closed  orbit  superim¬ 
pose  <1  on  the  surface  of  section  plot. 

Seven  separate  closed  curves  or  "loops"  are  also  clear¬ 
ly  evident  in  Figure  8  in  the  region  between  the  outermost 
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Figure  5.  Phobos  Surface  of  Section,  H  -6.832813 
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Figure  6.  P ho bos  Collision  Orbit 
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Figure  7.  Phobos  Surface  of  Section,  H  -6.8527 
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Figure  10.  Stable  Orbit  About  P ho  bos.  H  -h.8527 
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closed  curve  and  the  inner  closed  curves.  These  seven  loops 
are  the  chain  of  islands  associated  with  a  single  periodic 
orbit.  Seven  islands  exist  in  each  of  the  four  portions  of 
the  surface  of  section  plot.  The  orbit  traject..;.  associat¬ 
ed  with  this  resonance  travels  through  a  different  one  of 
the  seven  islands  for  each  encirclement  of  the  moon.  After 
seven  encirclements,  the  trajectory  comes  back  to  the  origi¬ 
nal  island  and  this  periodic  behavior  continues.  Figures  11 
and  12  show  the  periodic  orbit  superimposed  on  the  surface 
of  section  plot.  The  period  of  this  resonant  orbit  is 
approximately  seven  times  greater  than  the  period  of  the 
simple,  closed  periodic  orbit.  It  is  also  seven  times 
greater  than  the  orbit  period  of  Phobos  . 

Several  other  resonant  orbits  are  indicated  by  the 
surface  of  section  plots.  Figures  28  and  32  also  indicate 
the  presence  of  stable  periodic  orbits  with  a  resonance 
implied  by  seven  islands.  Figures  34,  38,  42,  44,  and  46 
indicate  even  greater  resonances.  Figures  34-36  represent 
an  especially  interesting  case  of  higher  order  resonance. 

Orbit  Evolution  .  As  the  Hamiltonian  is  increased, 
orbits  are  found  further  and  further  from  the  moon.  Trajec¬ 
tories  "on  the  edge"  of  the  stable  regions  in  the  surface  of 
section  plots  become  more  and  more  elliptical  and  their 
periapsis  points  move  closer  and  closer  toward  the  moon. 
The  orbits  become  less  and  less  stable  with  respect  to 
Phobos  and  are  essentially  orbits  of  Mars  instead,  which  are 
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perturbated  by  Phobos  . 

It  is  also  interesting  to  note  that  closed  periodic 
orbits  are  much  more  difficult  to  find  further  from  the 
moon.  The  closed  curves  associated  with  these  orbits  become 
very  narrow  and  elongated.  The  initial  conditions  for  these 
orbits  must,  therefore,  be  determined  very  precisely,  making 
them  much  less  navigable. 

The  characteristics  associated  with  the  evolution  of 
the  surface  of  section  plots  can  be  observed  in  Figures  13 
and  14.  The  upper  and  lower  branches  of  the  surface  of 
section  plot  in  Figure  13  clearly  demonstrate  that  the 
trajectories  are  almost  in  orbit  about  Mars.  The  surface  of 
section  points  on  these  branches  are  formed  as  the  moon 
rotates  beneath  the  very  elliptical  trajectories.  Figure  14 
implies  that  the  initial  conditions  for  the  stable,  closed 
periodic  orbit  indicated  must  be  determined  very  accurately; 
otherwise,  a  satellite  could  be  many  kilometers  from  its 
starting  point  after  completing  only  one  orbit  of  the  moon. 

Orbit  Direction.  The  surface  of  section  plots  provide 
no  indication  of  the  direction  of  the  orbits.  This  informa¬ 
tion  must  be  determined  from  other  a  priori  knowledge.  All 
of  the  stable  orbits  discovered  in  this  investigation  were 
retrograde  (counter-clockwise). 

Several  attempts  were  made  to  discover  prograde  (clock¬ 
wise)  orbits.  However,  the  test  trajectories  quickly  e  - 
caped  the  influence  of  the  moon  or  quickly  collided  with  the 


moon.  Few,  if  any,  useful  surface  of  section  points  could 
be  determined. 

Jefferies  discovered  in  his  work  that  stable  retrograde 
orbits  appear  for  much  smaller  values  of  the  Hamiltonian 
than  stable  prograde  orbits  (7:16).  The  Martian  moons  may 
be  too  small  to  support  any  prograde  orbits.  The  value  of  H 
for  which  stable  prograde  orbits  might  otherwise  appear,  may 
be  greater  than  that  required  for  escape  trajectories. 

De imos 

In  the  same  manner  as  the  Phobos  analysis,  several 
evenly-spaced  test  trajectories  were  initiated  along  the  X  - 
axis.  All  of  the  trajectories  were  integrated  for 
10,000,000  seconds  (116  days).  Because  the  orbit  periods 
were  generally  greater  than  those  of  the  Phobos  orbits,  a 
longer  integration  time  was  used.  The  longer  integration 
time  was  required  so  the  sufficient  number  of  surface  of 
section  points  needed  to  indicate  the  shape  and  size  of  any 
closed  curves  could  be  obtained. 

Table  II  provides  a  summary  of  the  characteristics  of 
each  of  the  stable,  closed  periodic  orbits  discovered  about 
Deimos.  The  orbit  period  of  the  first  orbit  listed  is 
approximately  one-fourth  of  Deimos's  orbit,  about  Mars.  As 
expected,  the  orbit  periods  increase  as  the  orbit  altitude 
is  increased. 

The  eccentricities  of  the  orbits  are  »pproximately  .33 
-  .35  with  the  exception  of  the  closer  orbits  which  are  more 


Table  II.  Characteristics  of  the  Stable  Orbits  About  Deimos 
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circular  in  nature. 


Table  II  also  indicates,  once  again. 


that  the  altitudes  of  the  orbits  increase  at  a  lower  and 
lower  rate  as  the  value  of  the  Hamiltonian  increases. 

The  initial  velocities  of  the  orbits  range  from  4.15 
m/s  for  the  lowest  orbit  (H  =  -2.738592)  to  43.8  m/s  for  the 
highest  orbit  (H  =  -2.7383). 

All  of  the  orbits  discovered  were  retrograde.  Attempts 
to  discover  stable  prograde  orbits  were  unsuccessful  as  in 
the  Phobos  case.  With  the  exception  of  the  low  altitude 
orbits,  the  evolution  of  the  surface  of  section  plots  was 
similar  to  that  seen  for  Phobos. 

Variations  in  the  semi-major  axis  that  are  quite  simi¬ 
lar  to  those  of  the  Phobos  orbits  were  discovered.  However, 
no  resonant  orbits  were  discovered.  This  may  result  from 
the  fact  that  Deimos  is  further  away  from  Mars  than  Phobos, 
thereby  limiting  the  resonances  associated  with  the  third 
body  effects.  At  low  altitudes,  orbits  representing  simple, 
rotating  Keplerian  ellipses  were  discovered.  None  were 
discovered  in  the  Phobos  analysis. 

Rotating  Ellipses.  The  surface  of  section  plots  for 
low  altitude  orbits  about  Deimos  are  quite  different  than 
those  of  the  other  orbits.  The  nearly  circular  closed 
curves  encompassing  the  entire  moon  in  the  surface  of  sec¬ 
tion  plot  of  Figure  15  indicate  the  presence  of  rotating 
orbits.  Figure  16  is  the  isolated  surface  of  section  formed 
by  only  one  trajectory  .  In  this  case,  the  orbit  is  a  Keple- 


Figure  15.  Deimos  Surface  of  Section,  H  -2.738592 
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rian  ellipse  that  rotates  about  the  moon.  The  periapses  and 
apoapses  occur  at  roughly  the  same  altitude  for  each  encir¬ 
clement  of  the  moon.  Therefore,  the  curves  traced  out  by 
these  points,  as  the  orbit  rotates  about  the  moon,  are 
nearly  circular.  Figure  17  demonstrates  this  behavior.  The 
orbit  superimposed  on  the  surface  of  section  of  Figure  16 
was  integrated  for  100,000  seconds. 

Figure  15  also  indicates  the  emerging  presence  of 
islands  associated  with  the  semi-major  axis  perturbation 
effect  discussed  earlier.  The  surface  of  section  for  this 
case  is  isolated  in  Figure  18.  The  stable,  closed  periodic 
orbit  associated  with  this  surface  of  section  is  displayed 
in  Figure  19.  With  the  exception  of  the  narrow  region  just 
discussed,  an  infinite  number  of  stable,  rotating,  periodic, 
elliptical  orbits  may  be  discovered  within  the  stable  region 
associated  with  the  closed  curves  of  the  surface  of  section 
plot  . 

The  discovery  of  rotating  ellipses  around  Deimos  and 
not  around  Phobos  is  most  likely  due  to  the  fact  that  Deimos 
is  much  further  away  from  Mars.  As  long  as  the  orbits  are 
very  close  to  Deimos,  the  moon's  gravitational  field  greatly 
dominates  the  possible  perturbation  effects  caused  by  the 
very  distant  planet. 

C_olli  s  i  on_  0 rjbi  ts_.  Trajectories  that,  collide  with 

Deimos  are  indicated  by  Figure  20.  One  such  trajectory  is 
shown  in  Figure  21.  All  of  the  periapsis  points  of  this 
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Figure  17.  Stable  Orbit.  About  Dp  i  mos  ,  H  r  -2.738592 
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rotating  orbit  are  below  the  moon's  surface. 

Orbit  Evolution.  The  evolution  of  orbits  about  Deimos 
as  the  Hamiltonian  and  altitude  increase  is  very  interest¬ 
ing.  The  orbits  discovered  just  above  the  surface  of  Deimos 
are  all  rotating  ellipses  (see  Figure  20).  However,  as  the 
altitude  of  the  orbits  increase,  variation  in  the  semi-major 
axis  begins  to  appear  (see  Figure  15).  The  areas  of  the 
surface  of  section  plots  Indicating  rotating  ellipses  de¬ 
crease  (see  Figure  22)  and  then  begin  to  disappear  complete¬ 
ly.  In  the  surface  of  section  of  Figure  23,  the  surface  of 
section  points  for  the  unstable  trajectory  initiated  just 
beyond  the  stable  region  were  included  to  indicate  the 
transition  to  non-rotating  orbits  with  varying  semi-major 
axes.  The  scattered  points  from  the  unstable  trajectory 
indicate  the  presence  of  a  disappearing  rotating  ellipse. 
When  H  is  further  increased,  only  non-rotating  orbits  are 
indicated  from  the  surface  of  section  (see  Figure  24).  The 
islands  first  appear  to  be  very  isolated.  Howe'-r,  as  the 
Hamiltonian  (and  orbit  altitude)  is  increased,  these  regions 
become  elongated  and  begin  to  approach  one  another  (see 


Figure 

25)  . 
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Cone  1  us i ons  and  Recommendations 

The  surface  of  section  technique  was  successfully  used 
to  determine  several  stable  periodic  orbits  about  both 
Phobos  and  Deimos.  Most  of  the  orbits  displayed  a  variation 
in  the  semi-major  axis.  Only  the  orbits  very  close  to 
Deimos  were  determined  to  be  simple  rotating  ellipses. 

Closed  periodic  orbits  are  much  more  difficult  to 
discover  further  away  from  the  moon.  In  addition,  the 
initial  conditions  associated  with  these  orbits  must  be 
determined  much  more  precisely  than  those  for  the  low- 
altitude  orbits.  The  gravitational  attraction  of  Mars 
dominates  for  the  high-altitude  orbits.  Therefore,  the  use 
of  orbits  in  close  proximity  to  the  moons  is  recommended 
over  the  distant  orbits.  Long  term  stability  may  be 
achieved  as  long  as  the  orbit  velocity  can  be  maintained  to 
within  a  few  meters  per  second. 

These  orbits  could  prove  useful  as  parking  orbits  for 
any  future  manned  or  unmanned  missions  to  Mars.  Because  of 
the  low  gravity  of  Phobos  and  Deimos,  the  landing  of  a 
vehicle  on  their  surfaces  from  a  low  altitude  orbit  would 
require  little  expenditure  of  fuel.  The  potential  extrac¬ 
tion  of  water  from  the  moons  makes  them  attractive  targets 
for  exploration. 

A  possible  extension  of  this  research  would  involve  the 
expansion  of  the  dynamics  model  to  include  lesser  order 
effects  such  as  the  n  o  n  spheric  a  1  characteristics  of  Mars 


and/or  the  effects  associated  with  other  bodies, 

the  sun . 


especially 


66 


Appendix  A:  Problem  Parameters 
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Appendix  B:  Equations  of  Motion 
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Appendix  C:  Phobos  Surface  of  Section  Plots 
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Figure  27.  Phobos  Surface  of  Section,  H  =  -6.8527 


71 


-120  -110  -100  -90  -80  -70  -60 

X  (km) 


Figure  30.  Phobos  Surface  of  Section,  H  =  -6.8526 


74 


400  -300  -200  -100 


II 


-J400  -300  -200  -100  0  100  200  300  400 

X  (km) 


Figure  31.  Phobos  Surface  of  Section,  H  =  -fe.8525 
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Figure  34.  Phobos  Surface  of  Section,  H 
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Figure  35.  Phobos  Surface  of  Section,  H  =  -6.8524 
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Figure  37.  Phobos  Surface  of  Section,  H  =  -6.8.r23 
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Figure  39.  Phobos  Surface  of  Section,  H  =  -6.8522 
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Figure  42.  Phobos  Surface  of  Section,  H  =  -6.8521 
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Figure  43.  Phobos  Surface  of  Section,  H  =  -6.852 
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Figure  46.  Phobos  Surface  of  Section,  H  =  -6.85 
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Figure  48. 
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53.  Deimos  Surface  of  Section,  H  =  -2.738589 
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Figure  55.  Deimos  Surface  of  Section,  H  =  -2.738587 
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